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Let v be a vector field in a bounded open set G C M d . Suppose 
that v is observed with a random noise at random points Xi , i = 
l,...,n, that are independent and uniformly distributed in G. The 
problem is to estimate the integral curve of the differential equation 

^-=v{x(t)), t>0,x(0)=x £G, 
at 

starting at a given point x(0) = xq £ G and to develop statistical 
tests for the hypothesis that the integral curve reaches a specified set 
L C G. We develop an estimation procedure based on a Nadaraya- 
Watson type kernel regression estimator, show the asymptotic nor- 
mality of the estimated integral curve and derive differential and 
integral equations for the mean and covariance function of the limit 
Gaussian process. This provides a method of tracking not only the 
integral curve, but also the covariance matrix of its estimate. We also 
study the asymptotic distribution of the squared minimal distance 
from the integral curve to a smooth enough surface FcG. Building 
upon this, we develop testing procedures for the hypothesis that the 
integral curve reaches F. 

The problems of this nature are of interest in diffusion tensor 
imaging, a brain imaging technique based on measuring the diffu- 
sion tensor at discrete locations in the cerebral white matter, where 
the diffusion of water molecules is typically anisotropic. The diffu- 
sion tensor data is used to estimate the dominant orientations of the 
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diffusion and to track white matter fibers from the initial location 
following these orientations. Our approach brings more rigorous sta- 
tistical tools to the analysis of this problem providing, in particular, 
hypothesis testing procedures that might be useful in the study of 
axonal connectivity of the white matter. 

1. Introduction. Let G C W 1 be a bounded open set. Suppose a vector 
field v : G \—> M. d is observed at points Xj G G, i = 1, . . . , n, with random errors; 
that is, the observations are 

V i = v(X i ) + Zi, 

where £,£i,--.,£n are i.i.d. bounded random vectors (r.v.) E£ = and 
Cov(£,£) = £. 

We are interested in the Cauchy problem for the differential equation 
(ODE) 

(1.1) ^ = v(x(t)), t>0,x(0)=x eG, 

at 

which of course can be equivalently written in an integral form, 

x(t) = xq + / v(x(s))ds. 
Jo 

Our goal is to provide an estimate X(t),t > 0, of its solution based on the 
data (Xi,Vi),i = l,...,n, and, most importantly, to study the asymptotic 
behavior as n — oo of statistics such as infn<t<T d 2 (X(t), T), where r C G is 
a given subset of G (most often, it will be the boundary of a specified region 
in G) and d(x, T) := inf{|x — y\ : y G T} is the usual Euclidean distance from 
i to T. This would allow us to suggest tests of the hypothesis that the true 
trajectory x(t), <t<T, reaches a certain region in G. 

Our main interest in this problem is related to its potential applications 
to diffusion tensor imaging (DTI), a technique in brain research introduced 
several years ago and often combined with conventional MRI (see, e.g., [7]). 
The diffusion of water molecules at a given location is characterized by a 
symmetric positive definite 3x3 diffusion matrix (diffusion tensor). The 
principal eigenvector of this matrix shows the dominant direction of the dif- 
fusion. In cerebral white matter, the diffusion is typically anisotropic and 
DTI allows one to recover its dominant directions by measuring the diffusion 
tensor field within voxels at a discrete set of locations and computing princi- 
ple eigenvectors of diffusion matrices (thus transforming the tensor field into 
a vector field, see Figure 1). The fiber tract then can be reconstructed by 
following the directions of the vectors in small steps from a specified initial 
location, which essentially means solving numerically the ODE generated by 
the vector field. This provides a noninvasive approach to study the axonal 
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Fig. 1. Shows the 3D vector field based on DTI data. The left graph is a fractional 
anisotropic (FA) map; the right graph gives the 3D visualization of the vector field inside 
a rectangular region of the FA map. 

connectivity of white matter fiber inside a brain region. The method is often 
referred to as white matter fiber tractography. 

Since the diffusion tensor field is being measured at a discrete set of lo- 
cations and each matrix in the field represents an average within a voxel 
corrupted with noise, it becomes crucial to use some methods of smooth- 
ing of tensor or vector fields or regularization techniques that restrict fibers 
to smooth paths. For instance, Basser et al. [2] applied B-spline smoothing 
to the tensor field; Poupon et al. [18] used Markov random field models to 
obtain a regularized estimate of the vector field. However, even fiber track 
estimates involving smoothing would possess a certain degree of variabil- 
ity and very little is known about quantitative ways to assess the variabil- 
ity of fiber track estimates (although Parker, Barker and Buckley [14] and 
Jones [10] suggest some Monte Carlo and bootstrap approaches to the prob- 
lem), which would facilitate the development of more rigorous approaches 
to statistical analysis of DTI data (the situation is somewhat different in 
conventional MRI and fMRI where approaches based on rather deep statis- 
tical understanding of the problem are becoming more common; see, e.g., 
[19] and [17]). This seems to be an important task especially because the 
mathematical models used in fiber tractography are rather involved and 
the existing methods utilize tools coming from very different areas (see 
[1, 2, 3, 6, 12, 13, 15, 16, 18, 20]). 

The goal of the paper is to make the first (and rather modest) step to- 
ward better theoretical understanding of statistical problems in DTI. Our 
approach to estimation of x(t),t > 0, is based on the Nadaraya— Watson type 
kernel regression estimate (NWE) V(x),x £ G, of the vector field v(x),x G G, 
which then is substituted for v into the ODE (1.1). The solution of the re- 
sulting Cauchy problem is an estimate X(t),t > 0. (This approach is some- 
what akin to that of Basser et al. [2] : the difference is that we are applying 
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smoothing to the vector field and not to the tensor field; also, we are using 
kernel regression based smoothing instead of splines.) We establish in Sec- 
tion 2 (the proofs are given in Section 3) the asymptotic normality of this 
estimator, that is, the weak convergence (in the space of continuous func- 
tions) of the properly normalized deviation process X(t) — x(t),t E [0,T], to 
a vector- valued Gaussian process on [0,T] with mean and (matrix-valued) 
covariance function that depend on the vector field v, on the covariance 
matrix E of the noise £ and on the kernel of NWE. We derive differential 
equations for the covariance function of the limit process, which allows us to 
develop a technique of simultaneous tracking of the fiber path and its covari- 
ance (see Section 4). We also study (Section 2) the asymptotic distribution 
(as n — > oo) of the distance info<t<T d 2 (X(t), T) from the estimated integral 
curve X(t),t £ [0,T], to a set F C G. The asymptotic distributions for such 
distances happen to be especially simple in the case when the minimum of 
the function [0, T] 3 t \— > d(x(t),T) is attained at a single point. In this case, 
the distributions are either normal or x 2 -type, and they depend on the ge- 
ometry of r and on whether the true integral curve x(t),t G [0, T] reaches 
r and in which way. These results allow one to bring into the analysis of 
DTI data some tools of rigorous statistical inference. In particular, one can 
use the asymptotic normality of X{t) to construct confidence ellipsoids for 
x(t) for a fixed t; one can go further and try to use the results on Gaussian 
processes to develop nonparametric confidence bands and hypotheses tests 
for the whole integral curve x(t),t £ [0,T]; one can develop statistical tests 
for the hypothesis that the true integral curve x(t),t S [0,T], reaches a spec- 
ified subregion of G; one can develop confidence intervals for the distance 
from x(t),t £ [0,T], to a subregion. The last two possibilities are especially 
important since they are related to the problem of axonal connectivity which 
is one of the central issues in DTI. In Section 4, some of these possibilities 
are studied both for simulated and for real data. 

There is a number of issues that (in our view) go beyond the scope of the 
paper, but that needs to be addressed to develop a comprehensive methodol- 
ogy based on our approach. First of all, the choice of NWE of the vector field 
v is relatively arbitrary. Similar theory could, in principle, be developed for a 
number of other smoothing techniques. Moreover, it might be more natural 
and statistically more appealing to do smoothing of the underlying tensor 
field and only then to compute the principal eigenvectors creating a vector 
field. However, methods of perturbation theory will be needed to develop the 
asymptotic theory of such estimators, which would make the mathematical 
analysis of the problem more involved. Also, it is not common in DTI (at 
least, to our knowledge) to measure direction vectors at random locations, 
so, regression models with fixed design would be more appropriate than the 
model with random design (which is used in the paper primarily because 
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the mathematics looks nicer). The vectors Vi are usually unit eigenvectors 
of diffusion matrices, so, it would be more natural to consider nonparamet- 
ric regression models for directional data rather than additive noise models. 
We are not exploring many important aspects of kernel-type nonparamet- 
ric estimation (such as data-driven choice of the bandwidth, minimax lower 
bounds, optimal convergence rates, adaptation, etc.). Finally, in fiber trac- 
tography it is of great importance to take into account the possibility of 
fiber paths branching or intersecting one another. This is not covered by our 
model (because of the uniqueness of the solution of ODE) and the extension 
of our results to this case poses some nontrivial problems. 

Realizing the importance of all these and some other issues, we, however, 
believe that the results we obtained so far might be of some interest for 
further development of a comprehensive statistical theory of DTI. 

2. A kernel estimate of integral curves and its asymptotic normality. 

We will assume that G C M d is a bounded open set of Lebesgue measure 1 
and, for simplicity, that Xi are i.i.d. uniformly distributed in G and that 
the r.v.'s {gi } are independent of {X^}. We also assume that supp(v) := 
{x : v(x) / 0} C G, which allows us to set v = outside of G. Furthermore, we 
need a smoothness assumption on the vector field v. Unless stated otherwise, 
we assume that it is twice continuously differentiable. 

We will use the following NWE of the vector field v (see, e.g., Efromovich 
W): 

i=i 

with a kernel K satisfying standard assumptions, in particular, 

/ K(x)dx = l, [ K(x)xdx = 0, 

and with some bandwidth parameter h = h n . It will be also convenient to 
assume that K has bounded support, where it is twice continuously differ- 
entiable (the last assumption can be replaced by a more mild one in most 
of the results, but it is not of great importance in the context of the paper). 
As a result, the estimate V{x) = outside a bounded neighborhood of G. 
Comparing with the standard NWE, our estimate is simpler: since the dis- 
tribution of Xi is known (it is uniform), we do not need to use the kernel 
density estimator in the denominator of V. In the case of nonuniform design 
some other smoothing techniques, such as local polynomial models (see [5]), 
might also be of interest. 

Then, we define a plug-in estimate of the solution x(t),t > 0, as the solu- 
tion X(t) = X n (t),t > 0, of the Cauchy problem 



(2.1) — j^ = V(X(t)), t>0, X(0) = x Q eG, 
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which is equivalent to the integral equation 

(2.2) X(t) = x + [ V(X(s))ds. 

Jo 

Note that since both v and V vanish outside a neighborhood of G (v 
actually vanishes outside G itself), x{t) and X(t) will remain in this neigh- 
borhood for all t > 0. 

To be specific, we assume that all vectors are vector columns; the sign * 
will denote transposition of vectors or matrices. Whenever it is convenient, 
we use the notation (•, •) for the inner product in W d and I denotes the iden- 
tity matrix. Also AA(/i, X) denotes the normal distribution with mean \x and 
covariance matrix X and Z ~ S) denotes a r.v. with this distribution. 

In what follows, we also need an estimate of the derivative of v and we 
use for this purpose 

1=1 

Under the above assumptions, V,V' are consistent estimates of v,v' uni- 
formly in R rf (see Lemma 1 below). 

Our first goal is to show that under the assumptions and nh d+s — > 

P > the sequence of stochastic processes -Jnh d - l (X(t) - x(t)),0 <t<T, 
converges weakly in the space C[0, T] = C([0, T], W 1 ) of Revalued continuous 
functions on [0,T] to the Gaussian process £(i) satisfying the SDE 

d£(t) = lL f K{u){v"{x{t))u,u)dudt + v'{x{t))£{t)dt 

(2.3) 2 hd 

+ Wv(x(t)))[Z + v(x(t))v*(x(t))]) 1/2 dW(t) 

with initial condition £(0) = 0, where W{t),t > 0, is a standard Brownian 
motion in M. d , 

ip(v):= / V(vT)dT, $(y) := / K(z)K(z + y)dz. 

Note also that v"(x(t)) involved in (2.3) is a d x d x <i-tensor and (v" (x(s))z, z) 
is a vector- valued quadratic form. In what follows, Mg(t) denotes the mean 
and C(t) the covariance matrix of (which does not depend on (3). In 
Section 4 we provide differential equations for Mp{t) and C(t). 

Theorem 1. Suppose that h n — ► and nh d+2 — > oo as n — ► oo. Then for 
allT>0 

sup \X n (t) — x(t)\ — ► asn— > oo, 

0<t<T 
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in probability. Suppose also that ra/i^ +3 — > (3 >0 as n — > oo. Let T > and 
suppose that for some 7 = jt > and for all < s <t <T 

1 



t 



t 

v(x(X))dX 



>r 



Then the sequence of stochastic processes \J nhn 1 (X n (i) — x(t)),0 <t<T, 
converges weakly in the space C[0,T] to the Gaussian process £(i),0 <t<T. 

Remark. The condition of boundedness of the support of the kernel K 
can be replaced by the conditions that the functions ^ and A (defined in the 
proof of Theorem 1) are integrable. The proof of Theorem 1 goes through 
and the theorem applies to such kernels as the Gaussian. 

We turn now to some consequences of the asymptotic normality. In par- 
ticular, we are interested in asymptotic properties of statistics of the type 
inf tg [ 0i T] Tn{d{X(t),r)), where T is a subset of G, d(x,T) is a distance from 
x to r and m is a monotone function [e.g., m(u) = u 2 or m(u) = u, u > 0]. In 
other words, we want to study the asymptotic behavior of the minimal dis- 
tance from the estimated integral curve X to a target set T. Such results are 
of statistical significance since they allow one to develop tests for hypothe- 
ses that the true integral curve x(t) is passing through a given region or to 
construct confidence intervals for the distance to the region. We will study 
this problem under the assumption that the function (p(x) := m(d(x,T)) is 
smooth enough, which leads to a somewhat more general question about 
convergence in distribution (subject to a proper normalization) of the se- 
quence inf t6[0jT] ip(X(t)) - inf t6[0jT] <p(x(t)). 

Theorem 2. Let x(t),t > 0, be an integral curve starting at x(0) = xo £ 
G. Suppose that (p:G*-^>~R is continuously differentiable. Denote 



M := jr € [0, T] : V (x(r)) = q mf y <p(x(t)) 

and suppose that M C (0, T). Finally, suppose the conditions of Theorem 1 
hold. Then the sequence of r.v.s 



nhr, 



inf u>(X(t))- inf <p(x(t)) 
te[o,T] te[o,T] 



converges in distribution to the r.v. inf Tg ^ ^(t)*</?'(x(t)). Ln particular, if 
the minimal set M consists of only one point r G (0, T), then the above 
sequence is asymptotically normal with mean Mg(r) and variance a 2 = 
(if 1 '(x(t)))* C '(t)ip' '(x(t)) . Suppose now that ip is twice continuously differ- 
entiable. If , for allreM, (p'(x(r))=0 and ip" (x(t))(v(x(t)) , v(x(t))) > 0, 
then the sequence of r.v.s 



nh n 



-1 



inf u>(X(t))- inf <p(x(t)) 
Ug[o,t] te[o,T] 
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converges in distribution to the r.v. 



2 T&M 



^'(x(r))«(r),f(r)) 



p"{x(t))(v(x(t)),v(x(t)))_ 
If the minimal set consists only of one point r, then the limit becomes 



ip"(x(r))(Z,Z) 



( V "(x(r))(v(x(r)),Z)f 



(P"(x(t))(v(x(t)),v(x(t))) _ 

Z ~N(M(3(t),C(t)) inR d . 

On the other hand, if for all u S M. d , ip" (x(t))(v(x(t)),u) = 0, then the dis- 
tributional limit of the sequence 



nh 



d-i 



inf (p(X(t))- inf ip(x(t)) 
te[o,T] te[o,T] 



is i inf rg Af ip" (x(t))(^(t), £(t)), which in the unique minimum case is ^"(x(r)) x 
(Z,Z). 

Consider two typical examples. 

Corollary 1. Let a £ G and x(t),t > 0, be an integral curve starting 
at x(0) = xq 6 G. Suppose that for some r £ (0, T) info<t<T — a| 2 = 
|x(r) — a| 2 , and, moreover, suppose thatr is the only point where the infimum 
is attained. Suppose also the conditions of Theorem 1 hold. If x(r) ^ a, then 
the sequence 



\fnht 1 



inf \X(t)-a\ 2 - inf \x(t) - a\ 

0<t<T 0<t<T 



is asymptotically normal with mean 2M / g(r)*(x(r) — a) and variance a 2 = 
A(x(t) — o)*C(t)(x(t) — a). Ifx{r) = a, then the sequence nh d ~ l 'mio<t<T \X (t) ■ 
a\ 2 converges in distribution to the r.v. 

(v(x(r)YZ) 2 



\z\ 



Z ~N{M p {t),C{t)) in 



\v(x(r))\ 2 ' 

Corollary 2. Let T := {x : \x — a\ = r} C G be a sphere. Let 
d(x, r) := inf \x — y\ = \ \x — a\ — r\ 

be the distance from x to V and let x(t),t > 0, be an integral curve starting 
at x(0) = xq £ G. Suppose that for some r G (0,T) 



inf d 2 (x(t),T) =d 2 (x(T),T) =:D 2 , 

0<t<T 
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and, moreover, suppose that r is the only point where the infimum is at- 
tained. Suppose also the conditions of Theorem 1 hold. If D 2 > 0, then the 
sequence 



nhn~ 



inf d 2 (X(t),T)- D< 

0<t<T 



is asymptotically normal with mean 2DMp(r)*n(x(T)) and variance 
a 2 = AD 2 n(x(T))*C(T)n(x(T)), 

where n(x) := ifEfr- If D 2 = and, moreover, the vector v(x(r)) is tangent to 

r, then the sequence nh^ -1 info<t<r d 2 (X(t), V) converges in distribution to 
the r.v. 7 2 , where 7 is a normal random variable with mean Mp{t)*ti{x{t)) 
and variance n{x{r))* C {t)ti{x{t)) . 



Remarks. 1. The result can be extended to more general smooth sur- 
faces r. In this case, n(x) would be the unit normal vector to T at the point 
x' G r that is the closest to x (assuming the uniqueness of such a point). 

2. Suppose H C G is an open nonempty subset of G with boundary 
dH = r. Let x(t),t G [0, T] , be the integral curve with initial condition x(0) = 
x ,x (£ H U T. If for some t £ [0,T] x{t) G H, then mi < t < T d 2 (x(t),T) = 
since x(t),t G [0,T], is a continuous function. Also, it easily follows from the 
first statement of Theorem 1 that with probability tending to 1 we have 
info<t<T d 2 (X(t), T) = (since X, being close to x uniformly in [0,T], must 
enter the set H and hence cross its boundary V). As a result, for any se- 
quence a n — > 00 



inf d 2 (X(t),T)- inf d 2 (x(t),r) 

0<t<T 0<t<T 
tends to in probability and in distribution. 



a n inf d 2 (X(t),T) 
0<t<T y w ' ' 



3. Proofs of the main results. We will need the following quite standard 
statement which is given without proof. 

Lemma 1. Suppose that h — > and nh d+2 — > 00 as n — > 00. Under the 
assumptions above, we have in probability 

sup \ V(x) -EV(x)\^0, 

x€R d 

sup |V^(2;) — v(x)\ — > and sup | V (x) — v'(x)\ — > 0. 
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Proof of Theorem 1. Consistency. Let y(t) :=X(t) — x(t). We have 

y(t)= f\v(X(s))-v(x(s))]ds 
Jo 

(V - v)(X{s)) ds + f\v(X{s)) - v(x(s))] ds, 
Jo 

which implies (using a Lipschitz condition on v and the fact that both X 
and x remain in a bounded neighborhood of G) that with some constant L 
for all te [0,7] 

\y(t)\ < T sup \V(x) - v(x)\ + L [ \y(s)\ ds. 

xGM d JO 

Using Lemma 1 and this Gronwall-Bellman inequality (see [8]), this easily 
implies consistency. 

Asymptotic representation of X — x. We will establish that 

(3.1) X(t)-x(t) = z(t) + 5(t), t£[0,T], 

where z(t) = z n (t),5(t) = 5 n (t) are sequences of stochastic processes such 
that the sequence converges in distribution to the 

Gaussian process £(t) and 

The following representation is obvious: 

y(t)= f'[V(X(s))-v(x(s))]ds 
< 3 ' 3 > Z 

(V -v)(x(s))ds+ [ v'(x(s))-y(s)ds + R(t), 



(3.2) sup \6(t)\=o p 

0<t<T 





where the remainder is defined as 



R(t) := f\v(X(s)) - V(x(s)) - v'(x(s)) ■ y(s)} ds 
Jo 



(V - v)(X(s)) -(V- v)(x(s))] ds 
- f\v(X(s)) - v(x(s)) - v'{x(s)) ■ y(s)} ds. 



Note that 

\(V-v)(X(s))-(V-v)(x(s))\ 
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(V - v)'(aX(s) + (1 - a)x(s)) da ■ y(s] 



n 



< sup \(V- v y(aX(s) + (l-a)x(s))\-\y(s)\ 

0<a<l 

< sup \V (x) — v (x)\\y(s)\. 



Also, 



\v(X(s)) - v(x(s)) - v'(x(s))y(s)\ 



[v'(aX(s) + (1 — a)x(s)) — v'(x(s))] da ■ y(s) 



< sup \v'(aX(s) + (1 — a)x(s)) — v'(x(s))\ ■ \y(s)\ 

0<a<l 

<r(\y(s)\).\y(s)\, 

where r(5) := sup xg]Rt i sup| y | <( 5 \v'(x + y) — v'(x)\ — > as 5 — > 0, since v' 
uniformly continuous on G. Then for all t S [0, T] 



(3.4) \R(t)\ < ( sup \V'(x) -v'(x)\ +r( sup \y(s)\ 

\0<s<T 



\y(s)\ds. 



Since sup < s <r \y(s)\ — > in probability and by Lemma 1 sup xgK d | V^ 7 (a?) 
v' (x) | — > in probability, we have for all T > 



(3.5) sup \R(t)\=oJ f T \y(s)\ds), 

0<t<T \J0 J 

which also implies 
(3.6) 



sup \R(t)\ =Op( sup \y(t)\ 

\0<t<T 



0<t<T 



Denote by z(t),t>0, the solution of the integral equation 

z(t):= f\v(x(s))-v(x(s))]ds+ [ v'(x(s))z(s)ds. 



In other words, z satisfies the equation 
dz{t) 



dt 



V(x(t)) - v(x(t)) + v'(x(t))z{t), z{0) = 0. 



Let S(t) := y(t) - z(t). We have 



S(t) = / v'(x(s))-S(s)ds + R(t), 
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which implies 

\5(t)\<\R(t)\+ /V(*(*))H*(*)I*»< SU P f \v'(x(s))\-\6(s)\ds. 

JO 0<t<T JO 

Applying again the Gronwall-Bellman inequality, we get 
\S(t)\< sup \R(t)\exp{ [ \v'(x(u))\du\ <C sup \R(t)\, < t < T, 

0<t<T [JO J 0<t<T 

with some constant C > 0, since the exponent is bounded. As a result, by 
(3.5), 

sup |<5(£)|=Op( / \y(t)\dt\ asn^oo. 
o<t<T \Jo ) 

Since y(t) = z{t) + 5(t), we also have 



sup \6(t)\ =o p ( [ \z(t)\dt 
o<t<T \Jo 



as n — > oo. 



It follows from the weak convergence of Vnh d 1 z n to the continuous stochas- 
tic process £ (established below) that 

(3 - 7) .f |2(t)l<ffi = >bi=T> 

Therefore, we also have (3.2). 

Consider now the following differential equation in M. d : 

t^i^m +Axmu{t) , „«,)= o, 

where F is a continuously differentiable function from [0,T] into M. d with 
F(0) = 0. This equation has the unique solution u(t) = UF(t), where U is 

a linear mapping from the space C^[0,T] of all continuously differentiable 
functions F on [0, T] with F(0) = into C[0, T]. Another routine application 
of the Gronwall-Bellman inequality shows that IA is a continuous (in fact, 
even Lipschitz) mapping with respect to the uniform distance. 
Denote 

rj n (t) := Vnhd- 1 [ (V(x(s)) -v(x(s)))ds, 



£ n (t) := Vnh^Znit), t>0. 

Then £ n =Uij n . We will show that the sequence of stochastic processes r/ n 
converges weakly in the space C[0,T] to the stochastic process r\ satisfying 
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the SDE 

dr)(t) = —— I K(u)(v" (x(t))u,u) dudt 

+ (iP(v(x(t)))[Z + v(x(t))v*(x(t))]) 1/2 dW(t) 

with initial condition r/(0) = 0. (Here and in what follows we often do not 
write the limits of integration for the integrals over M or over M. d .) Since 
U is a continuous mapping from (C^ [0, T], || • H^) into (C[0,T],|| • \\oo), 
this implies that £ n = lAr\ n also converges weakly to the Gaussian pro- 
cess Ur], which satisfies (2.3) and the initial condition £(0) = 0. Hence, 

Computation of mean and covariance. Define ft(s) := Im^s). Then 



rj n (t) = Vnh^ 1 J f t (s) ■ [V(x(s)) - v(x(s))] ds 

= Vnh^l J (x n {ft) ~ J Ms) ■ v(x(s)) ds^j , 



where 



X n (f) ■= I f(s)-V(x( S ))ds 

1 " r f x ( s )-X j 



_Lg| /(5)ir (4fti) rf „. W x j)+& ). 



3=1 

Let L denote the set of all bounded functions /onR such that the support 
of / is a subset of [0,T] and / is continuous almost everywhere in M. Note 
that L is a linear space and ft & L, t G [0, T]. In asymptotic representations 
for the expectation and covariance of X n (f), we assume that / G L. We start 
with KX n (f): 



EX n (f) = ± J fis^K ^ 8 ^ X y v (X)ds 

' f{s)jK(^pLy v{y)dyds 



f(s) J K(z) ■ v(x(s) — zh) dz ds 
f(s) J K(z) v{x(s))-hv'{x(s))-z 

+ ^ {v "( x (s)) z , z ) + o(h 2 ) 



dz ds 
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f(s)-v(x(s))ds-h J f(s)-v'(x(s)) ■ (^J K(z)zdz^j ds 

, 9 

+ I /(«)• / K(z){v"(x(s))z,z)dzds + o(h 2 ), 



2 

where we have used the substitution z = x ( s >~ y . Note that under the as- 



sumption nh d+3 — > f3 > we have 



ET?n(i)^— t [ K(z)(v"(x(s))z,z)dzds = Er](t). 



2 

Also, 



Cov(X n {f),X n (g)) = -l^Cov(| f{s)K( ^ {s) h ds ■ (v(x)+0, 



Since 



Gov (/ f(s)K i,. V (X),J 9{ s)K (5<4Z*) * . t ) = 0, 

0av(//(.)ir(5^)d..{,/»(.)Jr(2^)4-.m)-a 



we have 



Cov(X n (f),X n (g)) 
= (I) + (II) 

To handle (I), we write 

1 w/f^X^)-^^ £ ,* [„(x{u)-X 



1 

7^ 
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Note that 



E<^ K 



x(s)-X 



h 



K 



h 



h d / K(z)K[z + 



x(u) — x(s) 
~h 



dy 

dz 



x(u) — x(s) 



h 



Changing variable as u = s + rh, we get 



(3.9) (I) 



n 



h d - 



x(s + rh) — x(s) 



h 



f(s)T,g(s + rh) dr ds. 



Note that (x(s + rh) — x(s))/h — > v(x(s)) as n — > oo and also for all r and 
a.s. for s g(s + rh) — > g(s) as n — > oo (recall that the functions f,g £ L and 
hence are continuous a.e. in M). By assumptions K has bounded support, 
implying that the support of ^> is also bounded. At the same time, we have 



0<7< 



1 



u — s 

Therefore, the function 



v(x(X))d\ 



r i— >\I/(t)= sup VPIt 



0<s<u<T 



< sup \v(x)\ < +oo. 



v(x(X))d\ 



u — s 



also has bounded support and, since it is bounded, it is integrable. Thus, we 
can use Lebesgue dominated convergence to prove that 



x(s + rh) — x(s) 



h 



f(s)T,g(s + rh) dr ds 



*(v(x(s))T)dTf(s)Zg( S )ds= / ^(v(x(s)))f( S )Eg(s)ds, 



which along with (3.9) yields 

(i) 1+o(1) 



n 



h d ~ x 



ijj(v(x(s)))f(s)T l g(s)ds. 



[Indeed, the integration with respect to s is in a finite range, the function 
(s, r) i— ► f(s)Hg(s + rh) is uniformly bounded and 



x(s + rh) — x(s) 
h 



<*(r), 



16 V. KOLTCHINSKII, L. SAKHANENKO AND S. CAI 

so the dominated convergence can be used under the assumption that VP is 
integrable in R.] 

Similarly, the expression (II) can be written as 



x(u)-X 



1 

n 

Note that 



h 

f(s)-v(x(s))ds I g(u) ■ v (x(u)) du{\ + o(l)) 



v(X)-v*(X) \g(u)dsdu 



= h d J K(z)k(z + ^ - '^ juliM - zh)-v'(x{s) - zh)dz, 
where we use the substitution z — x ^"l~ y ,dy — h' l dz. Therefore, 
_L E (/ !(a)K (ffelz*) is . V{X) . V . (X) J K (&£*y u ) in 

= i// /<s) / k{z)k { z + iW r w )° (x(8) ■ zh) 

x v* (x(s) — zh) dz g(u) ds du, 
which after the change of variable u = s + rh becomes 

-^// /w /^( 2 + ^±I^W)„ w „)- 2ft) 

x v* (x(s) — zh) dzg(s + rh) dr ds. 

As before, we use the dominated convergence (under the same conditions) 
to show that the last expression is equal to 

l + o(l) 



nh d - 1 

= l + o(l) 
nh d ~ x 

implying that 



f(s) I K{z)K{z + Tv(x{s)))dzdTv{x{s)) ■ v* (x{s))g(s) ds 
f '(s)ip(v(x(s)))v(x(s)) ■ v*(x{s))g{s) ds, 



(II) = / 
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Finally, the covariance 

Cav(X n (f),XM) 
(3.10) = (I) + (II) 

= J i>{v{x{s)))f{s) • [£ + v(x(s)) ■ v*(x(s))] ■ g(s) ds. 



Thus, 



Cov(j7 n (£i),77 n (i 2 )) -> Cov(??(ti),r/(i 2 )) 



+ u(x(a)) • v*(x(s))] ds. 

Convergence of finite dimensional distributions (J.d.d.). We now turn to 
the proof of asymptotic normality of rj n (t),0 < t < T, in the sense of conver- 
gence of finite dimensional distributions. First we show that for all / G L 



VnhFi\X n (f) - J f(s)v(x(s)) ds 
converges to a normal distribution. Since by (3.8) 
VnhF^Ux n {f)- [ f(s)v(x(s))ds 



/(*)• / K(z){v"(x(s))z,z)dzds, 
it is enough to establish the CLT for 

Vn~h^(X n (f) ~ EXn(/)) = -r^=, E(Xi - E^-), 

vnh a+L - =1 

where 

xr-= J /(*)*( a(a) ~* J ' ) ds-(v(X 1 ) + t J ). 

Under the assumptions we have made it is easy to check Lyapunov's condi- 
tions for the CLT, and to this end we bound the fourth moment of Xj, 

|4 _ tcv, .*,. \2 



E\ Xj \* = E( x * jXj y 



x(s)-Xj\ ( x{si)-X 3 



2 

X (v(Xj)+S j y(v(X j )+£ j )f(8 1 )d8d8 1 ' 
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Under the assumption that v and £ are bounded, this gives with some con- 
stant C > 0, 



E| x /<Ce(/ //v(^^)/.(^-^)|/( s )||/( Sl )l^^ 

x \f(s)\\f(s 1 )\\f(s 2 )\\f(s 3 )\dsds 1 ds 2 ds 3 . 
By change of variable, we then get 

E\ Xj \*<Ch d jf e K(z)K (z + X(S1) ~ x{s) ) K [z + X{S2) ~ X{S) ) 

xi^ ^^ ) cte|/(-)||/(-i)||/(- a )| 
x |/(s3)| ds\ ds2 ds 3 

= CW*W K(z) K (z + T 1 x{s + Tlh) - x{s) ) 
Jm. 5 V T\h J 

x(s + r 2 h) — 

7 2 /i 

x(s + 73/t) — x(s) 

T 3 h 

X ^1/(8)11/(5 + 71^)11/^ + ^)1 
x l/( s + r 3^)l ^ s <^ r 2 dr 3 . 

Denote 

x(si) — x(s) ' 



x K \z + 7 2 
x iq z + 7 3 



A(7i,7 2 ,7 3 ) :=sup J K{z)K\z + t{- 



S\ — s 



x K[ z + r 2 xiS2) - x{s) )K(z + r 3 xiS3) - xisY ) dz, 



S2 — S J V S3 — S 

where the supremum is taken over all s, s±,s 2 , s 3 E [0, T). It follows from the 
conditions that the function A is integrable in M 3 . 
As a result, we get 

n Xj f<ch d+3 j j J A(7 1 ,7 2) 7 3 )(y \f( S )\u s y^j i/( s +7^)i 4 d s y /4 

X (/ \f{s + T 2 h)\Us^' 4 
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/ \f{s + r 3 h)\ dsj dTidT 2 dr 2 



Ch 



d+3 



\f(s)\Us 



A (n , r 2 , r 3 ) dr 2 dr 3 



It follows that with some constant C 



C 



1 A . _ F |4< Cn^+ 3 



implying Lyapunov's conditions for the CLT. This shows the asymptotic 
normality of 



Vnh d -i(X n (f)-EX n (f)) and Vnhf^[X n {f)- f(s)v(x(s)) ds 



for all f £ L. Hence, if /i, . . . , / m E L (which is a linear space), the CLT holds 
for any linear combination of /i, ■ • ■ ,/ m - Using the standard characteristic 
function argument, this shows that the joint distribution of (X n (fi), . . . , X n (f m )) 
is also asymptotically normal. Applying this to f = ft proves the conver- 
gence of finite dimensional distributions (f.d.d.) of the stochastic processes 
ijn{t),0 < t < T, to f.d.d. of the Gaussian process Tj(t), <t<T. 

Asymptotic equicontinuity. To prove the weak convergence of the sequence 
of processes r] n (t),0 <t<T, in the functional space C[0,T], it remains to 
check the asymptotic equicontinuity condition. Since 

r, n (t) = ^nh d -\X n {f t )-EX n {f t )) + ^nh d ~ l Ux n {f t ) -J f t (s)v(x{s))ds 

and the bias term v 7 nh d ~ l (EJ n (/() — / f t (s)v(x(s)) ds) tends to Erj(t) uni- 
formly in t G [0, T] due to (3.8), we have to consider only the process ( n (t) '■= 
V nh d ~ 1 (X n (ft) — EX n (/;)). To this end, we bound the fourth moment of 
X n (f)-EX n (f), 



E\X n {f)-EX n (f)\ e - 



(3.11) 



1 



1 



E 



n 



4 h 4d 



n(n — 1) 



(E|x-E X | 2 ) 2 + nE|x-E X | 



As before, 



E|x-Ex| 



<E|x| 
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E 



K 



r (x{s)-X\ K (x{ Sl )-X 



/(*) 



x {v{X)+i)*{v{X) + i)f{s 1 )dsds 1 



<c 

<Ch d 



-x 



K 



K{z)K[z + 



x( Sl )-X 
h 



x(si) — x(s) 

h 



x \f(s)\\f(s + Th)\dsdr 

*(T)\f(s)\\f( S + Th)\dsdT 
1/2 



\f(s)\\f(s 1 )\dsds 1 
dz\f(s)\\f(s 1 )\dsds 1 
dz 



< Ch d+l 

< Ch d+1 

< Ch d+1 



\f(s)\ Z ds 



*(r)dT / \f(s)\ 2 ds. 



\f(s + Th)\ z ds) dr 



1/2 



Plugging the bounds on K\x — E%| 2 and on E|x — Ex| 4 in (3.11) yields with 
a large enough constant C > 0, 

E(Vn/ l d - 1 |X n (/)-EX n ,(/)|) 4 



n 2 h 2d+2 



/l/WI'*)' + ^ 



l/(*)l 4 * 



We will apply it to / := — ft 2 with ii,i2 6 [0, T]. It easily follows that with 
some L>0, J \f tl (s) - f t2 (s)\ 2 ds < L\h - t 2 \ and / \f tl (s) - f h (s)\ 4 ds < 
L\t± — *2 1 - Therefore we have (with some C > 0) 



nuh)-ut2)\ A <c 



1 



l*i-*ar + ^3=rl*i-*al 



r? 



which gives E|C„(ti) -Cn(t 2 )| 4 < 2C\h -t 2 \ 2 for |t x -t 2 | < If now A n 

is a maximal nh d-i -separated subset of [0, T], then the standard Kolmogorov 
type of chaining argument shows that for all e > 

(3.12) limlimsuppj sup \<n(h) - ( n (t 2 )\ > e\ = 0. 

8^0 n^co Ui,t 2 eA n ,\ti-t 2 \<S > 

Let Ti n be a mapping from [0, T] into A n such that Vi £ [0, T] : \t — TT n t\ < 
l/nh d ~ l . Using the definition of ( n (t), we easily get (with some constant 



ESTIMATION OF INTEGRAL CURVES 21 

C>0) 

|Cn(ti)-Cn(*2)| <CT^Jnh d ~ l sup \V(x)-EV(x)\\tx-t 2 \. 

x£R d 

Therefore, 

sup \Ut) - U*nt)\ < CT r }— sup \V(x)-EV(x)\. 
te[o,T] Vn/i" 1 xeiR d 

Using Lemma 1, 

(3.13) sup |C„(t)-Cn(*r„t)|=o P (l). 

te[o,T] 

It immediately follows from (3.12) and (3.13) that 

limlimsuppj sup \( n (ti) - ( n (t 2 )\ > e\ = 0, 

5^0 n^oo Ui,i 2 e[0,T]|ti-t 2 |<5 J 

which is the asymptotic equicontinuity condition for the process Cn- d 
The proof of Theorem 2 can be found in [11]. 
4. Numerical implementation and examples. 

4.1. Remarks on numerical implementation. It is not hard to show that 
the mean of the limiting Gaussian process defined by (2.3) can be written 
as Mp(t) = \f]3M(t), where M satisfies the differential equation 

(4.1) ^1 = v '( x (t))M(t) + K(z)(v"(x(t))z,z)dz, M(0) = 0. 

The covariance matrix C(t) of £(t) does not depend on (3 and it satisfies the 
ODE 

(4.2) + v'(x{t))C{t) + C(t)v'(x(t))*, 
C(0) = 0. 

It is also easy to derive partial differential equations for the covariance func- 
tion of the process £(i), but they are not used in what follows. 

We will use a simple Euler type method to solve the ODEs numerically 
(obviously, more sophisticated numerical methods can also be useful here, 
with potential improvement of the results). Let 5 be a step size. Then the 
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following recurrent relationships provide an approximation of equations (4.1) 
and (4.2): 

(4.3) X k+1 := X k + V(X k )S with X = x 0) 

(4.4) C : = 0, 

C k+1 := C k + S[rP(V(X k ))(t + V(X k )V(X k )*) 
+ V'(X k )C k + C k V'{X k )% 

where X := n _1 £)?=i(Vj — V r (Xj))(Fj — V"(Xj))* is an estimate of the covari- 

ance matrix E of the noise £j. Consistency of the estimator S easily follows 
from Lemma 1. In practice, the noise £j is not necessarily homogeneous and it 
might make sense to use localized versions of the above estimate. Obviously, 
the recurrent relationships (4.3) and (4.4) can be solved simultaneously, so, 
in fact, our approach is based on simultaneous tracking of the "fiber path" 
and its covariance matrix. We are doing this for k = 1, . . . , N, N := 

It easily follows from the definition of the function tp (see Section 2) that if 
the kernel K is spherically symmetric, then ijj is also spherically symmetric, 
so it is a constant on the unit sphere in R . In applications, the vector field 
v consists of unit vectors. Hence, for a spherically symmetric kernel K, the 
^-factor in the ODE for C(t) and in (4.4) can be replaced by a constant, 
simplifying the equations. In what follows, we use the standard Gaussian 
kernel K, which of course is spherically symmetric. 

To estimate the function M(t), one needs an estimator of v". This can be 
done, for instance, by utilizing kernel estimators one more time. The entries 
of the estimator V" (which is a d x d x d-tensor) are defined as 

v" (x)- 1 T &2K ( x ~ Xi ) 
v M x )- nhd+2 X, dX]dxk { ~ h ) v i > 

Vp\ l = l,...,d, being the components of the vector V^. The kernel K can be 
taken to be the same as in the estimate of V, but the bandwidth parameter 
h = h n is different (so V" is not the second derivative of V). To make V" 
a consistent estimator of v" the assumptions h — > and nh d+ ^ — > oo are 
needed. The second assumption does not hold for the bandwidth h needed 
in Theorem 1. If K is the standard Gaussian kernel, then 

The following recurrence relation [i.e. to be solved simultaneously with 
(4.3) and (4.4)] provides a numerical approximation of equation (4.1): 

(4.5) M k+l = M k + 5[V'{X k )M k + \W{X k )) with M = 0. 



W( X ) := J K( Z )(V"( X)Z , z) dz = * £ ( * 
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Solving (4.3), (4.4) and (4.5) yields numerical approximations of X(t), 
C(i), and M(t),0 <t<T, that can be now used to compute mini<fe<jv d 2 (Xk, T), 
which is a numerical approximation of info<t<r d 2 (X(t), T), for a given set 
r, and also to compute other quantities needed for implementation of testing 
procedures. If the above minimum is attained at k and f := k6, then f can 
be used as an estimate of r for which the minimal distance from the true 
integral curve x(t),0 < t < T, to T is attained. If such a r is unique (as was 
assumed in Corollaries 1 and 2), then it is not hard to show consistency of 
f (under proper assumptions on 5). This provides an approximation of the 
limit distributions in Corollaries 1 and 2. 

The above considerations allow us to implement the testing procedures 
based on Corollaries 1 and 2. For instance, in the case of Corollary 1, the 
test statistic is approximated by 

(4.6) k:=nh d ~ l min \X k - a\ 2 . 

y ' l<k<N l 1 

Given a significance level a £ (0, 1), the hypothesis that the integral curve 
x(t),0 <t <T, passes through the point a (against the alternative that it 
does not) is rejected if A > A Q , where A Q is determined from the equation 
P{A > Aq,} = a. Here 

A:=\Z\ 2 - {V [ Xk )* Z) \ Z~M(VPM t ,C t ) inR d . 

i i \ V (x % )\* {VP fc ' k> 

We assume that h = (^) 1 ^ d+3 ^ with (3 > 0; we set (3 = if h = h n is such 
that re/4+ 3 0. 

Corollaries 1 and 2 can be also used to derive asymptotic approximations 
of the power of the test and to study how it depends on D (the minimal 
distance from the true integral curve to V). For instance, in the case of 
Corollary 1, the power can be approximated by the expression 

■{nh d ~ l )- l l 2 K - (nh d - l ) l / 2 D 2 - 2^DM(T)*n(x(r)y 



1_$ 



2D(n(x(T))*C(T)n(x(T))) 1 / 2 



where $ is the standard normal distribution function, D 2 : = info<j<T \ x(t) — 
a\ 2 and n(x) := ifEfr- Replacing M(r), C(r) and x(t) by their "estimates" 
leads to the following expression describing the dependence of the power on 
the true distance D: 

{nh d ~ l )~ l / 2 k a - (nh d ~ i yi 2 D 2 - 2JpDM?n(X [ ) 
(4.7) 1-$( V 



2D(n(X k )*C jc n(X j( )y/ 2 

We are not addressing in any detail an important problem of choosing 
the bandwidth parameter h. For a fixed t and h = (^) l ^ d+ ^ the asymptotic 



24 



V. KOLTCHINSKII, L. SAKHANENKO AND S. CAI 



formula for the mean squared error matrix of X is (see Theorem 1) 
E(X(t) -x(t))(X{t) -x(t))* 

« n-W + V[C(t)f3- 4/{d+3) + M{t)M(tri3 4 '^] 

(note that the convergence rate n -4 /^" 1-3 ) is optimal in a minimax sense 
provided that the vector field v is twice continuously differentiable; this 
can be shown, e.g., following the approach of [9], Theorem IV.5.1). This 
immediately implies the following formula for the mean integrated squared 
error: 

E [ \X(t)-x(t)\ 2 dt 
Jo 

T Tr(C(t))di/3- (d - 1)/(d+3) 



n -4/(rf+3) 



o 



+ [ T Tr(M(t)M(t)*)dt(3 4 ^ d+ V 
Jo 



10 

which can be easily minimized with respect to (3, and the minimal point 
j3 can be estimated based on the data (using the estimates of C and M). 
Since one might be interested in optimizing not the global deviation of X 
from x but rather the distance from x to a set T (as in Corollaries 1 and 2), 
an alternative is to use the asymptotics of these corollaries rather than the 
global result of Theorem 1. For instance, based on Corollary 1, the following 
asymptotic formula might be used: 



E 



inf \X(t)-a\ 2 - inf \x(t) - a\ 2 

0<t<T 0<t<T 

» n" 4 /( d+3 ) [4(x(r) - a)*C(T)(x(r) - a)/T (d ~ 1)/(d+3) 

+ A(M(t)*(x(t) - o)) 2 /3 4/(d+3) ], 

which again can be easily minimized with respect to f3, and the minimal 
point Pi can be estimated based on the data. One can also try to develop an 
approach based on maximizing the power of the hypothesis tests considered 
above. 



4.2. Several experiments with simulated and real data. We turn now to 
some of the results of our experiments with simulated and real data. First, 
we simulated a vector field with circular integral curves (Figure 2). It was 
observed at a finite number of random points uniformly distributed inside a 
rectangular domain in M. 2 with random noise. We used NWE to smooth the 
vector field and then computed an estimate of an integral curve starting at 
a given point by solving numerically the ODE generated by the smoothed 
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True Vector Field 



Not£¥ Vector Field 




■ 

- 

- 

• ■ 



... r.t... , 



-a 



NW smocHhed VekKily and BslimalBd IrsjBclDry 




EsJimaled Confidence interval by *tic!e covanarca matrix 




Fig. 2. The top left figure shows the true vector field whose x- and y -coordinates are 



given by formulas v x 



yj x 2 +y 2 ' 



("the circular field" ). The top right figure 



represents the noisy vector field obtained by adding to the true field independent copies 
of 0.5/7, Z ~ Af(0,I) in R 2 . The bottom left figure shows the results of smoothing of the 
noisy vector field and integral curve estimation using NWE. The total number of points 
in the rectangle n — 322. The NWE was computed with h = 0.85 and the step size used in 
the numerical solution of ODE was S — 0.02. Finally, the bottom right figure depicts 95 
percent confidence ellipses along the integral curve. 



field using Euler's method. Simultaneously with tracking the estimate of the 
integral curve, we have also tracked the covariance matrix of the estimate 
and used it to plot the 95 percent confidence ellipses along the integral curve. 
The results are shown in Figure 2 (see [11] for more experiments of similar 
nature). 

Next, by means of Monte Carlo simulation we studied the accuracy of 
normal approximation of the distribution of the distance from the estimated 
integral curve to a given point or to a given sphere (see Corollaries 1, 2). 
To this end, we simulated the random points and the noisy vector field as 
in Figure 2 and computed the estimated integral curve based on NW re- 
gression smoothing. We repeated these simulations independently N = 2000 
times and each time computed the square of the distance D 2 to the point 
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Fig. 3. Shows a circular true integral curve and locations of points and balls of interest 
in a Monte Carlo study of the distribution of distances (see Figures 4 &nd 5 below). 

(0,2) (labeled with + in Figure 3). The squared distance from this point to 
the true integral curve was D 2 = 1. Each time we also computed the esti- 
mate a 2 of the variance a 2 (see Corollary 1) and the standardized version 
of D 2 , given by the expression ^P^(l) 2 — D 2 ) (recall that d = 2 in our case). 
The histogram of the last variable is shown in the top part of Figure 4 in 
comparison with the standard normal curve. The bottom part of Figure 4 
shows the results of a similar simulation in the case of the distance from the 
estimated integral curve to a sphere (a circle in our case; see Corollary 2). 
There is deviation of the histograms from normality that is quite understand- 
able for a number of reasons: the fact that we ignored the bias Mp in the 
normal approximation; in the case when D 2 = 0, Corollary 1 suggests that 
the asymptotic distribution should be of x 2_ type rather than normal and 
because of this for small values of D 2 one can start seeing some deviations 
from normality for a finite sample; the variance a 2 needed in normalization 
was replaced by its estimate a 2 ; the numerical approximation we are using 
to compute the distance has a certain impact on the distribution; and, last 
but not least, the sample size n in our simulations is rather small for this 
type of CLT (n = 77). The Kolmogorov-Smirnov test clearly shows that 
these deviations from normality are very significant (p < 0.001). However, 
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Fig. 4. The top figure presents the histogram of standardized minimal squared distances 
from the estimated integral curves to the point x — (0, 2) obtained by the Monte Carlo 
simulations (N = 2000); the bottom figure shows the histogram of standardized minimal 
squared distances from the estimated integral curve to the ball with center x — (0, 2) and 
radius 0.1 obtained again by the Monte Carlo simulations (N — 2000,). The histograms are 
compared with the standard normal distribution. 

when n = 500 the p-value of the test becomes of the order 0.0542 and for 
larger values of n the deviations from normality are no longer statistically 
significant. 

Quite similarly, Figure 5 shows histograms of squared distances from the 
estimated integral curve to a specified point (the top figure) or to a specified 
circle (the bottom figure) in the case when the true integral curve passes 
through the point or is tangent to the circle. In this case, according to 
Corollaries 1 and 2, the asymptotic distribution of the squared distance 
should be of x 2 -type. 

Next we studied the power of testing the null hypothesis that the integral 
curve passes through a specified point of interest. The test is based on the 
second statement of Corollary 1. The test statistic is A given by (4.6). The 
top part of Figure 6 shows the true integral curve and also ten points of 
interest: one of them is on the curve (so that the null hypothesis is satisfied 
for this point) and nine other points represent alternatives. We estimated 
this integral curve based on n = 77 observations of a noisy vector field as in 
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Fig. 5. The top figure presents the histogram of the minimal squared distances from the 
estimated integral curve to the point x — (0,3) obtained by the Monte Carlo simulations 
(TV = 2000,) in comparison with a -type curve based on the theory. Note that now the 
point is on the true integral curve. The bottom figure shows the histogram of the minimal 
squared distances from the estimated integral curve to the ball with center x — (0, 2.9) and 
radius 0.1 obtained by the Monte Carlo simulations (W = 2000J again in comparison with 
a x 2 -type curve based on the theory. The empirical distributions in this case are much 
closer to y^-type than the distributions shown in Figure 4- 



Figure 2. We repeated the experiment 1000 times, each time simulating the 
data, estimating the integral curve and testing the hypothesis with signifi- 
cance level a = 0.05. The red curve shown in the bottom part of Figure 6 
represents the empirical estimation of the power of our test (the frequency 
of rejecting the null hypothesis) for each of the alternatives. The blue curve 
represents the value of the power based on the theoretical formula (4.7) 
(which seems to consistently overestimate the power). 

Figures 7(a)-(c) give some examples of fiber tracking and visualization of 
95% confidence ellipsoids for real DTI data and Figure 8 represents what 
we call the p-value map. It can be used to assess the degree of connectivity 
of points with a given path. 

Although it is not the goal of this paper to develop a comprehensive 
methodology of statistical analysis of DTI data that takes into account 
more complicated issues, such as crossings or branchings of fibers, some 
of our results suggest possible ways to address these problems (as well as 
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Fig. 6. The top part shows the true integral curve and selected points of interest for 
measuring the power. The bottom part represents graphs of the power function based on 
Monte Carlo study (the red curve) and based on the theoretical formula (4-7) (the blue 
curve). 

the problem of stopping criteria of the tracking). Usually the vector fields 
involved in DTI problems satisfy the condition \v(x)\ = 1 since this is the 
field of unit principal eigenvectors of diffusion matrices. If the vector field is 
smooth around a given point x G 67, then V(x), as a weighted local average 
of unit vectors v (Xj) of approximately the same direction plus a small noise, 
should have norm close to 1 (for large enough n). On the other hand, if two 
or more trajectories intersect at the point x, then V{x) becomes a weighted 
local average of unit vectors of several different directions and, as a result, 
jV^x)! will be, with high probability, significantly smaller than 1. Based on 
these heuristics, one can try to test the null hypothesis that a trajectory 
x(t) starting at a given point xq does not intersect another trajectory at 
a given moment t > 0. It is natural to use the deviation jV^A^t))! 2 — 1 as 
a test statistic and to reject the null hypothesis if the value of this statis- 
tic is small enough. To make it more precise, one can use the following 
asymptotic result that easily follows from the CLT via a version of the delta 
method (and the result of Theorem 1, is also needed to handle the remain- 
der terms): under the assumptions of Theorem 1, if h n — > 0, nh^ +2 — ► oo 

and n/i^ +3 — > as n— > oo, then the sequence of r.v.s drilled V(X(t))\ 2 — 1) 
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Fig. 7. Illustrates the use of the proposed tracking procedure on real DTI data. (&) 
presents a single estimated fiber trajectory. The blue point shows the starting seed point. 

(b) shows the visualization of the 3-D confidence ellipsoid (C.E.) of the tracking procedure. 

(c) is an enlarged visualization of a 3-D fiber tracking trajectory in a given region of the 



converges in distribution to the normal random variable with mean and 
variance a 2 := 4/ K 2 (u) du(l + v(x(t))*T,v (x(i))), which can be replaced by 
a plug-in estimate a 2 in a straightforward way. The statistic 
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Fig. 9. Shows the results of tracking a branching trajectory (for simulated 2-D data). 
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becomes asymptotically standard normal and the test can be designed in 
the usual way. However, it should be pointed out that the convergence rate 
in this asymptotic result is rather slow and the impact of the remainder 
terms is significant. Moreover, when instead of the additive noise model 
one considers a more realistic model in which \Vi\ = 1 (this is the case in 
DTI applications), the asymptotic distribution of v n is no longer standard 
normal. 

The computation of the test statistic can be included in the fiber tracking 
algorithm. If the test detects an intersection at some moment t, the com- 
putation of the weighted local average V(X(t)) can be replaced by local 
clustering of vectors Vi around this point. For instance, one can use for this 
purpose a weighted version of k- means clustering. As a result, there will 
be more than one continuation of the trajectory from the point X(t). A 
simulated example of testing and tracking a branching trajectory using this 
approach is given in Figure 9. The branching point of the true trajectory is 
(0,0). The figure also shows the values of the test statistic at several points 
along the trajectory. 

A detailed discussion of this and other applications of our methodology to 
DTI and its comparison with other methods goes beyond the scope of this 
paper and will be given in further publications in more specialized journals 
on neuroimaging. 
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